%%% Written by  Daniel Lewis, 2020, for "Robust Inference in Models 
%%% Identified via Heteroskedasticity", Review of Economics and Statistics.
%
% This file conducts the simulation study underlying Figures 3-5 in the
% paper. MC_size.m should be run first to produce quantiles necessary for
% SAP calculations.

clear; clc
warning('off','all');
options=optimset('MaxFunEvals',10000,'MaxIter',10000,'TolFun', 1e-10,'Display', 'off','TolX',1e-10);
%% Setup

% Setting up parameters
T_emp=833;
Tc_emp=759;
params=[.7035 -.3098 .0039 0.0001 0.0070 0.0004]; % parameter values from empirical application (specification 1)
factor=params(5)/params(3); % proportional variance change for non-policy shock
dtilde_T=params(6)-factor*params(4); % equivalent to sigma^2_{2,c}*sigma^2_{1,p}/sigma^2_{1,c}*d/sqrt(T)
dtilde=dtilde_T*sqrt(T_emp); % equivalent to sigma^2_{2,c}*sigma^2_{1,p}/sigma^2_{1,c}*d - since we care about sigma^2_{2,p} as a function of d, not d itself, this is all we need
Tfrac=Tc_emp/T_emp; % fraction of control days in empirical sample
T=800; % sample size for simulation
Tc=ceil(Tfrac*T); % control sample size for simulation
Tp=T-Tc; % policy sample size
ident=[1/10,1,10]*dtilde; % degrees of identification for simulation
N=1000; % Number of Monte Carlo samples
H=[1,params(2);params(1),1]; % Empirical H matrix
V1=diag(params(3:4)); % Structural variances for control regime
alt=-3:.1:3; % grid to construct alternatives for H_12
axis=alt+params(2); % grid of alternatives for H_12
load('quanttab.mat'); % load quantiles for SAP calculation (produced by MC_size.m)

% Storage
Fstatpw=NaN(N,3,length(alt));
Sstatpw=NaN(N,3,length(alt));
eststore=NaN(N,3,length(alt));

% Compute symbolic Jacobian
Hs=sym('H',2); % symbolic H matrix
assume(Hs,'real') ;
Hs(1)=1;Hs(4)=1; % unit diagonal normalization
Vc=diag(sym('Vc', [2 1])); % control regime structural variances
Vp=diag(sym('Vp', [2 1])); % policy regime structural variances
moms=[reshape(Hs*Vc*Hs',4,1);reshape(Hs*Vp*Hs',4,1)]; % compute moments
moms([3 7])=[]; % eliminate redundant moments
theta=[Hs(2); Hs(3); diag(Vc); diag(Vp)]; % stacked parameter vector
D=jacobian(moms,theta); % compute symbolic Jacobian
Dt=D; % rescaling rows of symbolic Jacobian by sample share
Dt(1:3,:)=D(1:3,:)*Tc/T;
Dt(4:6,:)=D(4:6,:)*Tp/T;
grad=matlabFunction(Dt,'Vars',theta); % convert symbolic Jacobian to function for fast evaluation


%% Runs

for j=1:3
    
    fprintf('Identification strength: %d\n',j);
    drun=ident(j);
    V2=[params(5),0;0,factor*params(4)+drun/sqrt(T)]; % scaling variance by sqrt(T) according with local-to-unity device
    truen=[params(1:5),V2(end)]'; % set true parameter values under the null for this calibration - note this depends on drun
    
    for v=1:length(alt)
        fprintf('Alternative: %.2f\n',axis(v));
        true=[params(1) axis(v) params(3:5) V2(end)]'; % set true parameter values, same as truen above but with the alternative for H_12
        Htrue=[1,true(2);true(1),1]; % true H, incorporating alternative
        
        parfor i=1:N
            rng(i) % seed rng
            % storage
            etasq1=NaN(Tc,3);
            etasq2=NaN(Tp,3);
            
            % draw data
            eps1=mvnrnd([0;0],V1,Tc)'; % shocks in control regime
            eps2=mvnrnd([0;0],V2,T-Tc)'; % shocks in policy regime
            eta1=Htrue*eps1; % compute innovations
            eta2=Htrue*eps2; % compute innovations
            eta1=eta1-mean(eta1,2); % demean so outer product equals covariance
            eta2=eta2-mean(eta2,2); % demean so outer product equals covariance
            
            % estimate structural parameters
            for s=1:length(eta1) % compute outer product of innovations in control regime
                mat=eta1(:,s)*eta1(:,s)';
                etasq1(s,:)=[mat(1:2),mat(4)]; % store vech
            end
            for s=1:length(eta2) % compute outer product of innovations in policy regime
                mat=eta2(:,s)*eta2(:,s)';
                etasq2(s,:)=[mat(1:2),mat(4)]; % store vech
            end
            cov1=eta1*eta1'/Tc; % covariance in each regime
            cov2=eta2*eta2'/Tp;
            mat=cov2*cov1^-1; % ratio of covariance matrices
            [a,b]=eig(mat);  % estimate columns of H using Brunnermeier et al (2019) approach
            if b(1)>b(4) % if the first shock has the largest proportional variance change
                Hint=a(1,1)/a(2); % parameter of interest is 1,1 element, but need to re-order and impose unit-diagonal normalization
                Ho=a(2,2)/a(3); % doing the same for other free element in H
            else
                Hint=a(1,2)/a(4); % parameter of interest is 1,2 element, need to impose unit-diagonal normalization
                Ho=a(2,1)/a(1); % doing the same for other free element in H
            end
            Hhat=[1,Hint;Ho,1]; % estimated H matrix, correctly ordered and normalized
            n1=Hhat^-1*cov1*Hhat'^-1; % structural variance matrix in control regime
            n2=Hhat^-1*cov2*Hhat'^-1; % structural variance matrix in policy regime
            est=[Ho,Hint,n1(1),n1(4),n2(1),n2(4)]; % putting together estimated parameter vector
            
            % compute test statistics
            [~,~,omega]=momfast(etasq1,etasq2,T,est); % compute moment covariance matrix
            ghat=grad(Ho,Hint,n1(1),n1(4),n2(1),n2(4)); % evaluate Jacobian at estimated parameter values
            V=(ghat'*omega^-1*ghat)^-1; % asymptotic variance
            Fstatpw(i,j,v)=T*(est(2)-truen(2))'*V(2,2)^-1*(est(2)-truen(2))/1; % compute and store F-statistic for H_12
            try % first use true values as start values
                [est1,S1,flag1]=fmincon(@(p) momfastsub(etasq1,etasq2,T,p,H(1,2),2),truen([1 3:end]),[],[],[],[],[-Inf; eps*ones(4,1)],[],@ordering,options); % minimize objective function conditional on null value of H_12
            catch
                S1=NaN; est1=NaN; flag1=NaN;
            end
            try % also try unrestricted estimates as start values
                [est2,S2,flag2]=fmincon(@(p) momfastsub(etasq1,etasq2,T,p,H(1,2),2),est([1 3:end])',[],[],[],[],[-Inf; eps*ones(4,1)],[],@ordering,options); % minimize objective function conditional on null value of H_12
            catch
                S2=NaN; est2=NaN; flag2=NaN;
            end
            if flag1>0 && flag2>0 % if both find a valid minima, record the smaller
                Sstatpw(i,j,v)=min(S1,S2);
                if S1<=S2 % and save the corresponding estimate
                    eststore(i,j,v)=est1(1);
                else
                    eststore(i,j,v)=est2(1);
                end
            elseif flag1>0 % if only the first finds a valid minima, save its output
                Sstatpw(i,j,v)=S1;
                eststore(i,j,v)=est1(1);
            elseif flag2>0 % if only the second finds a valid minima, save its output
                Sstatpw(i,j,v)=S2;
                eststore(i,j,v)=est2(1);
            end
            
            if 100*floor(i/100)==i % report progress every 1000 samples
                fprintf('%d draws complete\n',i);
            end
        end
    end
end

%% Figures 3 and 4: Power curves
% Note that Figure 3 is merely a truncated version of Figure 4; this code
% produces Figure 4 in its entirety

figure
lines={'-','--','-.',':'}; % specify line styles
tit={'$\hat{d}/10$','$\hat{d}$','10$\times\hat{d}$'}; % titles
subplot1(3,1,'FontS',12,'Gap',[.01,.07]);

for j=1:3
    subplot1(j)
    hold on
    Fp=squeeze(Fstatpw(:,j,:)); % select slice of data
    Sp=squeeze(Sstatpw(:,j,:)); % select slice of data
    datamat=[nanmean(Fp>finv(.95,1,inf),1);nanmean(Sp>chi2inv(.95,1),1);nanmean(Sp>chi2inv(.95,6),1)]; % compute rejection rates relative to standard critical values
    for a=1:size(datamat,1)
        plot(axis,datamat(a,:)','LineStyle',char(lines(a)),'LineWidth',2)
    end
    ax = gca;
    set(gca,'TickLabelInterpreter','LaTex')
    set(gca,'FontSize',16)
    if j==3
        ax.XTick = [-3 -2 -1 -.3098 0 1 2];
        ax.XTickLabel = {'-3','-2', '-1', '$H_{12}$', '0', '1','2'};
    end
    xlim([min(axis) max(axis)]);
    xl=xline(-.3098); % add a line at null for H_12
    xl.LineStyle='--';
    xl.Color='k';
    grid on
    str=char(tit(j)); % add title
    title(str,'Interpreter','LaTex');
    hold off
end
subplot1(1) % Add legend to first plot
legend('Wald','K_{sub}','K_{proj}','Location','NorthWest','FontSize',12)
legend('boxoff')

%% Figure 5: SAP curves

figure
lines={'-','--','-.',':'}; % specify line styles
tit={'$\hat{d}/10$','$\hat{d}$','10$\times\hat{d}$'}; % titles
subplot1(3,1,'FontS',12,'Gap',[.01,.07]);

for j=1:3
    subplot1(j)
    hold on
    Fp=squeeze(Fstatpw(:,j,:)); % select slice of data
    Sp=squeeze(Sstatpw(:,j,:)); % select slice of data
    datamat=[nanmean(Fp>quanttab(j,3));nanmean(Sp>quanttab(j,4))]; % compute rejection rates relative to quantiles from size simulations (SA crit. values)
    for a=1:size(datamat,1)
        plot(axis,datamat(a,:)','LineStyle',char(lines(a)),'LineWidth',2)
    end
    ax = gca;
    set(gca,'TickLabelInterpreter','LaTex')
    set(gca,'FontSize',16)
    if j==3
        ax.XTick = [-3 -2 -1 -.3098 0 1 2];
        ax.XTickLabel = {'-3','-2', '-1', '$H_{12}$', '0', '1','2'};
    end
    xlim([min(axis) max(axis)]);
    xl=xline(-.3098); % add a line at null for H_12
    xl.LineStyle='--';
    xl.Color='k';
    grid on
    str=char(tit(j)); % add title
    title(str,'Interpreter','LaTex');
    hold off
end
subplot1(1) % Add legend to first plot
legend('Wald','K','Location','NorthWest','FontSize',12)
legend('boxoff')

%% Shock ordering constraint
function [c,ceq] = ordering(thetahat) % imposes ordering constraint: largest proportional variance change is in policy shock, epsilon 2
c = thetahat(4)/thetahat(2)-thetahat(5)/thetahat(3);
ceq = [];
end
